#!/usr/bin/env python

from math import *
import hooloovookit.astro.convert as convert
import hooloovookit.astro.geometry as geometry
import hooloovookit.astro.Sun as Sun
import hooloovookit.astro.Moon as Moon
import hooloovookit.astro.astrodatetime as astrodatetime

APOLAT = 32.7797556
APOLONG = 105.8198305

def lunskybright(alpha,rho,altmoon,alt):
    """ From Skycalc: Evaluates predicted LUNAR part of sky brightness, in 
        V magnitudes per square arcsecond, following K. Krisciunas
        and B. E. Schaeffer (1991) PASP 103, 1033.
        
        alpha = separation of sun and moon as seen from earth, in Degrees
        rho = separation of moon and object, in Degrees 
        altmoon = altitude of moon above horizon, in Degrees
        alt = altitude of object above horizon, in Degrees
        
        
        The original C code has the following extra parameters, taken here to be constants:
        kzen = zenith extinction coefficient
        moondist = distance to moon, in earth radii
        
        all are in decimal degrees. """
    
    if altmoon < 0.0: return 0.0
    
    kzen = 0.19 # Zenith extinction
    moondist = 60.27 # Earth radii
    
    rho_rad = radians(rho)
    alpha = 180. - alpha
    Zmoon = pi/2. - radians(altmoon)
    Z = pi/2. - radians(alt)
    moondist = moondist / (60.27)  # divide by mean distance
    
    istar = -0.4 * (3.84 + 0.026*fabs(alpha) + 4.0e-9*alpha**4.) # eqn 20
    istar =  (10.**istar) / moondist**2
    
    if fabs(alpha) < 7.: # crude accounting for opposition effect 
        istar = istar * (1.35 - 0.05 * fabs(istar))
    
    # 35 per cent brighter at full, effect tapering linearly to 
    #   zero at 7 degrees away from full. mentioned peripherally in 
    #   Krisciunas and Scheafer, p. 1035.
    fofrho = 229087. * (1.06 + cos(rho_rad)**2.)
    
    if fabs(rho) > 10.:
        fofrho = fofrho + 10.**(6.15 - rho/40.) # eqn 21
    elif (fabs(rho) > 0.25):
        fofrho = fofrho + 6.2e7 / rho**2 # eqn 19
    else:
        fofrho = fofrho + 9.9e8 # for 1/4 degree -- radius of moon! 
    
    Xzm = sqrt(1.0 - 0.96*sin(Zmoon)**2)
    
    if(Xzm != 0.):
        Xzm = 1./Xzm
    else:
        Xzm = 10000.
    
    Xo = sqrt(1.0 - 0.96*sin(Z)**2)
    
    if(Xo != 0.):
        Xo = 1./Xo
    else:
        Xo = 10000.
    
    Bmoon = fofrho * istar * (10.**(-0.4*kzen*Xzm)) * (1. - 10.**(-0.4*kzen*Xo))  # nanoLamberts
    
    if(Bmoon > 0.001):
        return 22.50 - 1.08574 * log(Bmoon/34.08)  # V mag per sq arcs-eqn 1
    else:
        return 99.

def mjdRADec2skyBright(mjd, ra, dec):
    dtObj = convert.mjd2datetime(mjd)
    
    moon = Moon.Moon(dtObj)
    moonRA, moonDec = moon.ra, moon.dec
    
    sun = Sun.Sun(dtObj)
    sunRA, sunDec = sun.ra, sun.dec
    
    # alpha
    moonSunAngle = geometry.subtends(sunRA, sunDec, moonRA, moonDec, units="DEGREES")
    
    # rho
    moonObjectAngle = geometry.subtends(moonRA, moonDec, ra, dec, units="DEGREES")
    
    moonAlt, moonAz = convert.raDec2AltAz(moonRA, moonDec, APOLAT, APOLONG, dtObj)
    objAlt, objAz = convert.raDec2AltAz(ra, dec, APOLAT, APOLONG, dtObj)
    
    if moonAlt > 0 and objAlt > 0:
        bright = lunskybright(moonSunAngle.deg, moonObjectAngle.deg, moonAlt, objAlt)
    else:
        bright = 0
    
    return bright

"""
if __name__ == '__main__':
    import sys
    
    if sys.argv[1] == "-h" or sys.argv[1] == "--help":
        print "This script computes the lunar sky brightness in V band magnitudes on a given MJD at a given RA,DEC."
        print "Usage: moonSkyMag <MJD> <RA> <Dec>"
        sys.exit(0)
    
    try:
        mjd = float(sys.argv[1])
        ra = float(sys.argv[2])
        dec = float(sys.argv[3])
    except IndexError:
        print "Usage: moonSkyMag <MJD> <RA> <Dec>"
        sys.exit(1)
    
    bright = mjdRADec2skyBright(mjd, ra, dec)
    
    # Compare with moonSkyMag.c
    print "brightness in V band: %.2f" % bright
else:
    raise ImportError("This module is an executable! Not meant to be imported.")
""" 

if __name__ == '__main__':
    import sys, argparse
    
    parser = argparse.ArgumentParser(description="Compute the lunar sky brightness in V band magnitudes on a given date at a given position on the sky.")
    parser.add_argument('datePos', metavar='', type=str, nargs='+',\
                        help='Date and position, in the format: <DATE> <COORD1> <COORD2>')
    parser.add_argument('-d', '--date-units', dest="dateUnits", default="MJD")
    parser.add_argument('-c', '--coordinate-units', dest="coordUnits", default="degrees")
    parser.add_argument('-e', '--equitorial', dest="equitorial", action="store_true", default=True)
    parser.add_argument('-g', '--galactic', dest="galactic", action="store_true", default=False)
    args = parser.parse_args()
    
    if len(args.datePos) != 3: raise ValueError("You must supply a date and two coordinates! e.g. subtends DATE COORD1 COORD2")

    try:
        x,y = float(args.datePos[1]), float(args.datePos[2])
    except:
        raise ValueError("Error parsing coordinates!")
        
    if args.coordUnits.lower() == "degrees":
        pass
    elif args.coordUnits.lower() == "radians": 
        x,y = degrees(x), degrees(y)
    else:
        raise ValueError("Coordinate Units specified, %s, are not valid! Use 'degrees' or 'radians'" % args.coordUnits)
    
    if args.equitorial:
        ra,dec = x,y
    if args.galactic:
        ra,dec = convert.galactic2RaDec(x, y)
    
    if args.dateUnits.lower() == "mjd":
        mjd = float(args.datePos[0])
    elif args.dateUnits.lower() == "jd":
        mjd = float(args.datePos[0]) - 2400000.5
    elif args.dateUnits.lower() == "iso":
        import dateutil.parser as dup
        dt = astrodatetime.datetime.fromDatetime(dup.parse(args.datePos[0]))
        mjd = dt.mjd
    else:
        raise ValueError("Coordinate Units specified, %s, are not valid! Use 'degrees' or 'radians'" % args.coordUnits)
    
    bright = mjdRADec2skyBright(mjd, ra, dec)
    
    # Compare with moonSkyMag.c
    print "brightness in V band: %.2f" % bright
    
else:
    raise ImportError("This module is an executable! Not meant to be imported.")